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Abstract. In this paper we investigate the effects of stochasticity in the spatial and temporal 
distribution of supernova remnants on the anisotropy of cosmic rays observed at Earth. The 
calculations are carried out for different choices of the diffusion coefficient D(E) experienced 
by cosmic rays during propagation in the Galaxy. The propagation and spallation of nuclei 
(with charge 1 < Z < 26) are taken into account. At high energies (E > 1 TeV) we assume 
that D(E) oc (E/ZY , with (5 = 1/3 and S = 0.6 being the reference scenarios. The large scale 
distribution of supernova remnants in the Galaxy is modeled following the distribution of 
pulsars with and without accounting for the spiral structure of the Galaxy. Our calculations 
allow us to determine the contribution to anisotropy resulting from both the large scale 
distribution of SNRs in the Galaxy and the random distribution of the nearest remnants. The 
naive expectation that the anisotropy amplitude scales as 5 a oc D(E) is shown to be a wild 
oversimplification of reality which does not reflect in the predicted anisotropy for any realistic 
distribution of the sources. The fluctuations in the anisotropy pattern are dominated by 
nearby sources, so that predicting or explaining the observed anisotropy amplitude and phase 
becomes close to impossible. Nevertheless, the results of our calculations, when compared to 
the data, allow us to draw interesting conclusions in terms of the propagation scenario to be 
preferred both in terms of the energy dependence of the diffusion coefficient and of the size of 
the halo. We find that the very weak energy dependence of the anisotropy amplitude below 
10 5 GeV, as observed by numerous experiments, as well as the rise at higher energies, can 
best be explained if the diffusion coefficient is D(E) oc E 1 /^. Faster diffusion, for instance 
with <5 = 0.6, leads in general to an exceedingly large anisotropy amplitude. The spiral 
structure introduces interesting trends in the energy dependence of the anisotropy pattern, 
which qualitatively reflect the trend seen in the data. The inhomogeneous spatial distribution 
of the sources in the Galactic disc induces a large scale anisotropy which is not sensitive 
to the stochastic nature of nearby SNRs: we find that this additional contribution to 5a 
becomes more important for large values of the size of the halo, H. The two terms are 
comparable in size for H ~ 2 kpc which corresponds to the scale height of the gradient of 
the spatial distribution of SNRs in the Galaxy. The dependence on energy of 5a{E) is close 
to monotonic when the large-scale, regular term dominates, and does not seem to reflect the 
observed anisotropy amplitude. Both contributions to the total anisotropy are illustrated 
and discussed with the help of semi-analytical results. 
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1 Introduction 

The anisotropy of cosmic rays (CRs) may be a precious tool to probe the propagation of 
CRs throughout the Galaxy. A comprehensive understanding of the origin of CRs can only 
be achieved through a combined investigation of the spectrum, chemical composition and 
anisotropy. In paper I [1] we have discussed the spectrum and chemical composition of CRs 
as observed at Earth when the sources are assumed to be supernova remnants (SNRs). The 
stochastic distribution of SNRs in the Galaxy was found to lead to interesting effects: if 
diffusion is parametrized in terms of a diffusion coefficient D(E) oc E s , the best fit to the all- 
particle spectrum is found to require 7 b s = j + 5 = 2.67, where 7 is the slope of the injection 
spectrum. In a rigidity dependent scenario for the CR acceleration and propagation, the 
knee in the all-particle spectrum is determined by the combination of the spectra of different 
chemicals. Qualitatively the resulting spectrum well reproduces the observed one, but the 
stochastic nature of the supernova events leads to fluctuations which imply that the naive 
expectation 7^ = 7 + 5 does not hold strictly, a fact that makes it less straightforward to 
infer the actual injection spectrum required to fit the data. In an approximate sense, however, 
it remains true that a stronger energy dependence of the diffusion coefficient corresponds to 
harder injection spectra. This is a very important point in that the best known acceleration 
mechanism, diffusive shock acceleration (DSA) [2-4], invariably leads to spectra very close to 
power laws with slope 7 ~ 2 (naively leading to require 5 ~ 0.6—0.7). In fact, modern versions 
of DSA, including non-linear effects induced by the pressure of CRs in the shock region 
and by the super- Alfvenic streaming of CRs in the upstream region make two fundamental 
predictions: 1) concave spectra, possibly even harder than E~ 2 at high energies, and 2) 
large values of the maximum achievable energy, plausibly as large as the knee for protons 
[5], and Z times larger for nuclei with charge Z. The two aspects are strictly connected with 
each other: efficient acceleration leads to a high CR pressure, which in turn is responsible 
for a pronounced precursor and concave spectra; large CR acceleration efficiencies also imply 
magnetic field amplification and possibly higher maximum energies if CRs manage to interact 
resonantly with the self-generated magnetic perturbations. 

It has recently been proposed that the non-linear theory of DSA (hereafter NLDSA) 
may lead to somewhat steeper spectra if the velocity of the scattering centers in the amplified 
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magnetic field is taken into account in solving the transport equation for CRs at SNR shocks 
[6-8]. In this case the spectra of accelerated particles may be less concave, and have an 
approximate power law shape with slope ~ 2.1 — 2.2. The steeper slope of the injection 
spectrum follows from both the finite velocity of the scattering waves and from the convolution 
over time of the acceleration history of SNRs throughout their evolution in the interstellar 
medium (ISM). However, there are two disappointing elements in this picture: first, these 
values of 7 would still imply 5 ~ 0.5 — 0.6, in order to ensure j b s ~ 2.7; in addition, the 
essential role played by the velocity of the scattering centers introduces much uncertainty 
since the results depend rather strongly on poorly understood details such as the wave helicity, 
as well as the reflection and transmission of waves at the shock surface (in principle the effect 
discussed here could even lead to harder spectra). Another disappointing ingredient, that 
has important consequences on the shape of the spectrum, is represented by the details of 
the magnetic field amplification, as recently discussed by [9]. This again goes in the direction 
of steepening the spectra, although a non linear analysis of this effect has not been carried 
out as yet. 

As shown in Paper I, one element that can help us discriminating the case of fast diffusion 
(5 ~ 0.6, 7 ~ 2.1) from slow diffusion (5 ~ 1/3, 7 ~ 2.3) is represented by the shape of the all- 
particle spectrum and the chemical composition: fast diffusion leads to sizeable fluctuations 
in the spectra, which woud imply that fitting the all particle spectrum depends upon the 
specific realization of the source distribution, invalidating any given physical criterion like 
7+5 = jobs- The case of fast diffusion also leads in general to negligible spallation of nuclei, so 
that the Helium spectrum at Earth is roughly parallel to that of protons, in contradiction with 
recent CREAM [10] and PAMELA [11] data that suggest a harder He spectrum above TeV 
energies. There might be different explanations for the hardening of the He spectrum (see e.g. 
[12, 13]), but it is however noticeable that if the diffusion coefficient is properly normalized 
and taken to scale with energy as E 1 / 3 , this puzzling feature may be very naturally explained 
as a result of propagation from discrete sources. 

In Paper I we also discussed how the CR spectrum at Earth is affected by the details of 
the spatial and temporal distribution of SNRs. We showed that taking into account the spiral 
structure of the Galaxy leads to small changes compared to the case when the disk is taken 
as smooth. The changes become increasingly appreciable the smaller the assumed size of the 
magnetized galactic halo. We also considered different scenarios in terms of escape of the CRs 
from the acceleration sites and the possibility that supernovae of type II may preferentially 
be clustered within OB associations/superbubbles. We found that the resulting all-particle 
spectrum is very little sensitive to these assumptions. 

While the changes induced in the spectrum by the choice of 5 or by the stochasticity in 
the spatial and temporal distribution of the sources are small, they become very large when 
the anisotropy amplitude is considered, a fact that makes this quantity very interesting to 
investigate. This is the purpose of the present paper. 

The two benchmark cases, 5=1/3 and 5 = 0.6 lead to very different predictions for the 
anisotropy, as already noticed in Ref. [14]: the anisotropy obtained in the case of fast diffusion 
(5 = 0.6) is in general much larger than the observed one. Moreover the fluctuations related to 
the specific realization of sources are very large, making it impossible to obtain a quantitative 
prediction to compare with data. The comparison must necessarily be qualitative. 

The fluctuations appear, though with different strengths, in all scenarios that we dis- 
cuss in the following, and they play a crucial role in explaining the observed anisotropy. The 
naive prediction of the diffusion model with a continuous distribution of sources is that the 
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anisotropy should increase with energy exactly as E s . This is the type of anisotropy signal 
that can be calculated in numerical approaches such as GALPROP or DRAGON, which 
however can only single out the anisotropy signal due to the large scale inhomogeneous dis- 
tribution of the sources (and/or to the off-center position of the solar system in the Galaxy). 
In fact, the anisotropy pattern of individual realizations shows a rather irregular behavior 
with local dips and peaks in some energy ranges, or sometimes showing approximate con- 
stancy over large energy intervals, similar to what the data show. The theoretically predicted 
scaling of the anisotropy with the diffusion coefficient can only be recovered when the size of 
the halo is taken to be larger than 4 kpc, in agreement with the fact that when the halo is 
large enough the large scale inhomogeneity of the source distribution becomes increasingly 
dominant over the small scale fluctuations due to nearby sources. 

The paper is organized as follows: in § 2 we summarize the Green function formalism 
introduced in Paper I; in § 3 we derive simple estimates of the effect of fluctuations on the 
spectrum for a homogeneous but discrete distribution of sources in the disc. Our results on 
anisotropy for different models of source distribution and escape mechanisms of CRs from 
the sources are presented in § 4. In § 5 we discuss the dependence of the expected anisotropy 
signal on the size of the halo. We conclude in § 6. 



2 Summary of the Green functions formalism 

In this section we briefly summarize the Green function formalism that has been introduced 
in Paper I. The diffusive transport of cosmic rays from a point source located at a position 
r s = (x s , Us, Zs) and injecting a spectrum N(E) at a time t s can be written as: 

dn k (E,r,t) _ „ rn fr , ,m r,sp fTP , lw -f .v . AT , mxu . U 3 W * 



Of 



V [D k (E)Vn k (E,f,t)} - T s k p (E)n k (E,f,t) + N k (E)5(t - t s )5 d (f- f s ), (2.1) 



where n k (E,f,t) is the density of particles of type k (nuclei) with energy E at the location 
r and time t, D k (E) is the diffusion coefficient assumed to be spatially constant and T sp (E) 
is the rate of spallation of nuclei of type k to lead to lighter nuclei. Each source is assumed 
to produce nuclei of H (jfe = 1), He (k = 2), CNO (k = 3), Mg-Al-Si (k = 4) and Fe (k = 5). 
All injection spectra are assumed to be in the form: 



N k (E) oc£- 7 exp 



E 



(2.2) 



where the normalization and the maximum energy of each species are chosen as discussed in 
Paper I. The maximum energies scale with the charge Z of the nucleus. 
The diffusion coefficient is taken in the form: 



D(E) = W™D 28 (^ v y C m 2 s- 1 (2.3) 

for rigidity R > 3 GV, the secondary to primary ratios and the abundances of unstable 
isotopes can be best fit by choosing D2s/H kpc = 1.33 for 5 = 1/3 and D<2,&/H kpc = 0.55 for 
5 = 0.6, where H kpc is the height of the halo in units of kpc. 

The diffusion region is assumed to be a cylinder of infinite radius and half height H, 
inside which the diffusion coefficient is spatially constant. The fact that the radius is taken 
potentially infinite is not a limitation to the calculation because our source distributions are 
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concentrated in a region of a few kpc radius (the Galactic disc) and for practical purposes 
the size of the diffusion region can be taken to be much larger than the radius of the disc. 
The escape of CRs from this model Galaxy occurs through the upper and lower boundaries. 
The escape is modeled by assuming that n k vanishes at z = ±H, so that the escape flux 
through the surfaces z = ±H is described by D k (E)-^-\ z= ±H- 

In this paper, as in Paper I, we concentrate our attention on particles with energy above 
~ 1 TeV, therefore in Eq. 2.1 we ignored terms of advection and terms of possible second order 
reacceleration which are both potentially important at much lower energies. We also neglect, 
in computing the flux of nuclei, the contribution of secondaries produced in spallation events. 
From the point of view of how diffusive propagation is treated, the calculation described here 
is not much different from GALPROP or similar propagation codes. The only two relevant 
differences are in the more phenomenological way in which spallation is described and in the 
fact that here we do not impose the free escape boundary condition at some finite radius along 
the lateral sides of the cylinder. The latter, as discussed above, is not a relevant limitation, 
although it might lead to some corrections for large values of if. As to spallation, its rate 
is defined in terms of the process cross section and of the density of gas in the diffusion 
region as explained in Paper I. This treatment is adequate as long as only primary nuclei are 
considered, as we do here, while a more refined treatment would be required if we were to 
describe rare nuclei, such as B and 10 Be or similar secondary products. 

The Green function that satisfies the correct boundary condition at z = ±H can be 
obtained through the image charge method and can be written as follows: 



Gk(r,t;f s ,t s ) = Nk ^ E ' t 1 S ), 2 exp [-T^(E)t] exp 

[A-KD k T\ ' 



(x - x s ) 2 + (y -y s 
4£>*t 



-l) n exp 



<) 5 



4A-T 



(2.4) 



where z' n = (—l) n z s + 2nH are the z coordinates of the image sources. It is easy to check that 
Qk{x, y, z = ±if, t; x s ,y s , z s ,t s ) = by simply expanding the sum term. In what follows, the 
Earth will be located at (x, y, z) = (R@,Q, 0), where Rq = 8.5 kpc is the distance of the Sun 
from the center of the Galaxy. It is however useful to keep Eq. 2.4 in its most general form 
because anisotropics are related to the spatial derivatives of the Green function calculated 
at the detection location. 

It is important to realize that the Green function formalism outlined here allows us to 
take into account a completely arbitrary spatial distribution of the sources in the Galaxy. 
The same is true for the temporal evolution of the injection of CRs by each source: one 
can introduce an arbitrary time dependence of the particle spectrum injected per unit time, 
Q(E,t), and simply integrate Eq. 2.4 over time. 

We already discussed the importance of this latter point in Paper I, where we worked 
out the expression of N(E, t) for the two scenarios that are most relevant for CR acceleration 
by SNRs. These are: 1) "burst injection", with all particles from a SNR injected at the same 
time with a spectrum N{E)\ 2) continuous injection, with Q(E,t) oc 5{E — E max (t)), namely 
with particles at the maximum energy escaping the SNR at all times after the beginning of 
the Sedov- Taylor phase, at time T$t- In this latter scenario, the maximum energy the SNR 
can provide is Em and is reached at the time T$t- Later on, the maximum energy decreases 
with time allowing particles of progressively lower energy to leave the system. The whole 
process lasts ~ (3 — 30) x 10 4 years, and the convolution over time of these peaked spectra 
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leads to a power law injection spectrum [15] which is not directly related to the spectrum of 
particles at the shock. The global injected spectrum from an individual SNR is likely to be 
the result of the superposition of the spectrum of particles escaping the SNR from upstream 
and of those that escape after the end of the Sedov- Taylor phase. The latter component, 
however, is heavily affected by adiabatic energy losses, and therefore it does not contribute 
to the energy region around the knee. 

In the first scenario, namely for "burst injection", the Green function in Eq. 2.4 is 
already the solution that we are seeking, with 



N k (E) 



ilk ~ 2)r]k£ki 



E o,k 



1 



-7/C+2 



E \ " 7fc 
E-o) 6XP 



E 



E 



max,k 



(2.5) 



where rj k is the fraction of the kinetic energy of the blast wave, e k i n , that goes into accelerated 
particles of type k. The reference energy Eq^ is assumed to be 1 GeV for protons (for 
heavier nuclei its numerical value is not really important since we simply rescale the injected 
protons spectrum in order to fit the spectra observed at Earth). In Eq. 2.5, E maXyk is the 
maximum energy of particles of type k. In the expression above and in what follows we 
always assume that the injection spectrum is steeper than E~ 2 , since flatter spectra would 
result in unreasonable choices for the energy dependence of the diffusion coefficient (leading 
to exceedingly large anisotropy and even to breaking the regime of diffusive propagation, see 
below) . 

For the second model of injection, we assume that the maximum energy of the acceler- 
ated particles increases during the ejecta dominated phase of the remnant evolution to reach 
a maximum at the beginning of the Sedov- Taylor phase, and start decreasing afterwards. 
Particles at the instantaneous maximum energy escape the SNR at all times t > T$t [15], 
carrying a fraction rj(t) of the explosion energy. We take the decrease with time of E max 
in the form of a power-law, E max (t) = Ej[,f(t/TsT)~ a , with a > chosen in such a way as 
to guarantee that E max (jsNR) — Eq, where tsnr is the time at which the SNR dies out 
and we take Eq = lGeV for protons. With these assumptions, and using the Sedov- Taylor 
description of the shock dynamics we find for the injected spectrum: 



Q(E,t) = V (t) 



1 



Em Tst 



Tst 



a-l 



t>T ST 



(2.6) 



We further assume r](t) = ??o(*/^sr)^ , with j3 > 0, as discussed in Paper I. Supplied with 
this latter condition, the expression in Eq. 2.6 can be used in Eq. 2.4. Integration over time 
then leads to (see Paper I for details): 



n(E,f,t) = rjo e ki 



6\° _L_ 

5/ aE 



M 



[4itD(E)t*} 3/2 



exp [-Tf(E)r*] exp 

+oo 

E ("l) n exp 



Em \ ( Em 
E ) \ E 

(x - x s ) 2 + (y- y s ) 
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21 



(2.7) 
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where 



T* = t-t a - T ST 



E J 



From Eq. 2.7 it is clear that the solution is non-zero only for 

i 

Em 



1 < 



E 



< Min 



1 + 



T"SNR t — t s 



Tst Tst 



(21 



Of all the inequalities in Eq. 2.8, the interesting case is that of recent SNRs, for which 
(t — t s )/TsT < 1 + tsnr/Tst- In this case, at time t, the SNR only contributes particles 
with energies 

't-U 



E > Em 



Tst 



(2.9) 



This lower limit on energy comes from the fact that lower energy particles are still confined in 
the source. This feature is peculiar to the scenario of continuous injection and is completely 
unrelated to propagation effects, while the latter set, in all cases, an additional, independent 
bound on the minimum energy of CRs reaching us from a given source. 



3 Anisotropy and its fluctuations 

The main goal of the present paper is to discuss the implications of diffusive CR propagation 
from discrete sources on the anisotropy signal measured at Earth. In a regime of diffusive 
propagation, the anisotropy in a given direction is defined as 

5 _ 3D{E) V s n CR {E,r,t) 
c n C R 

where ncn(E) is the CR number density as measured at the position of the Earth. The 
density gradient, as calculated over many realizations of sources, has two contributions, one 
that is sensitive to the mean distribution of sources over large spatial scales, and one that 
is sensitive to fluctuations induced by random nearby sources. Whether one or the other 
dominates the anisotropy signal depends in general upon the spatial scales on which the 
gradients in the source distribution appear. 

It is easy to see that for a homogeneous distribution of sources in the disc, the mean 
value of Sg vanishes in all directions, simply because there are no gradients in the disc 
(z = 0) in such a situation, except for boundary effects close to the edge of the disc. Even 
in this situation, however, the measured anisotropy will be non-zero, due to the small scale 
inhomogeneity of the source distribution in the solar neighborhood. It is clear that this 
contribution will not show the scaling with D(E) that one expects based on Eq. 3.1. For 
an inhomogeneous distribution of sources, as the one that we think appropriate to describe 
SNRs in the Galaxy, it is not straightforward to predict what the anisotropy will be as a 
function of energy, since this depends on whether the small or large scale contribution will 
dominate. In order to quantify these effects we follow here the argument of [16] , already used 
in Paper I to calculate the fluctuations in the spectral shape of CRs at Earth as due to the 
stochastic distribution of sources. 

In order to write the anisotropy it is useful to introduce the particle current. In the 
case of diffusive propagation, this reads: 

Jcdr, t) = -DVn CR {r, t) . (3.2) 
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The overall anisotropy will then result from the sum of two terms, representing the large-scale 
and small-scale contribution respectively. These can be written in terms of Jcr as: 

ftu = (3-3) 
±cn C R{E) 

and 

, (S Jcr -S Jcr) 1 / 2 

0A2 = J ~=T • (3.4) 

^cn C R{E) 

If we consider now M independent sources, each producing Nq particles, following the 
formalism of [16] and of Paper I, we can write the average particle current entering Eq. 3.3 
as: 

(Jc R (r,t)) = N AfD(E) J d 3 ^dt'V r g(r,t-y,t')P(^,t') , (3.5) 
while for its fluctuations, entering Eq. 3.4, we have: 

(SJ OR (r, t)SJ CR (f, t)) = N 2 MD(E) 2 T WV r Q(f, t; r", t") ■ V r , Q(f ', t'; f", t")P(f" ', t") . 

Jd 

(3.6) 

In the previous expressions P(r, t) is the probability of having a source at the position 
(r,t) at time t, Q(r,t;f*,t') is the Green function for transport of particles from (f,t') to 
(f,t) and terms of 0(1/ M) have been neglected. 

Now let us consider a homogeneous spatial distribution of the sources in the disc, de- 
scribed by the probability function 

P(r,t) = ^2- (3-7) 
nK d 

In this case, it is easy to show that: (JcR,(r, t)) = and hence the large scale term in Eq. 3.3 
is zero. As to the small scale term, for purpose of illustration, let us compute it for the case 
of protons alone, so that spallation can be ignored, and let us consider an infinitely thin 
Galactic disc, so that (5 Jcr ■ SJcr) = {8JcR,x8JcR,x) + (5JcR,y&JcR,y)- This case is clearly 
much simpler than the scenarios we actually considered in our calculations, but it still allows 
us to discuss a few key features. From Eq. 3.6 we have: 



K f°° dr f Rd r 2 



(SJcrSJcr) = -m I tjziz^ I j2Drf exp 

{l 2 + m 2 )H 



vRd 2 J Tmm (4vrdr)3 J {2DtY 

r n2 i ™2\rr2 



r 



2 



2Dt 



Dt 



(31 



x E E(-!)' +mex p 

m=— oo l=— oo 

where both integrals above can be calculated analytically, and, if we assume H 2 / D(E)T m i n S> 
1, the only important term in the sums is easily seen to be that corresponding to I + m = 0. 
The result for the rms current is then: 

1 1ZN 2 

(SJcrSJcr) - m ^£^2 ■ (3-9) 
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At the same time, it is straightforward to calculate the average CR density (see Paper I). 
For H <C Rd, this can be written as: 



N(E)K H 

{nCR) = 2,R\ D(E) ■ (3 - 10) 

Using the expressions in Eq. 3.9 and Eq. 3.10, from Eq. 3.1 we finally obtain: 

3 D{E)^ (R d \ 1 
Sm ~ 2*W CUV* {~Hj T~- (3 - U) 

First thing to notice in Eq. 3.11 is the strong dependence of the small scale anisotropy 
(in terms of normalization and in terms of energy dependence) on the cutoff time T m i n , which 
we were forced to impose in order to avoid the divergence in the integral over r. One can 
see that, if T m i n is naively chosen as a given number, the anisotropy is a very slow function 
of energy, 5a2 ~ E & / 2 . On the other hand, the time T m j n should carry information about 
the closest, most recent sources around the observer. For instance the time T m i n could be 
interpreted as the time over which one source goes off within a distance from Earth such that 
CRs from that source reach us within a time T m i n . This condition is expressed as 



-P- [4irD(E)T min ] =l-> T min 



AKD(E) 



R 2 



-1/2 

(3.12) 



With this prescription the small-scale anisotropy reads [14]: 



r 3 1 D(E) 

= W^"^' (3 - 13) 

which shows an energy dependence Sa2 ~ E s . Moreover Sa2 D(E)/H, namely a scaling 
that is the inverse of that of the CR flux, ncn(E) cx H/D{E). The fact that both quantities 
depend upon the ratio of D and H, and that the same ratio at, say, 10 GeV, is fixed by 
measurements of the B/C ratio, implies that the assumed size of the halo has little effect in 
changing the small-scale anisotropy or the spectrum. However this statement depends on the 
specific choice of T m j n , that we have made and cannot be taken as general: when fluctuations 
are important, the actual anisotropy can easily be expected to depend strongly upon the 
specific distribution of nearby and most recent sources. 

On the other hand, for any realistic distribution of sources in space, the zero-th moment 
of the anisotropy, Sai, does not vanish and corresponds to a large scale anisotropy due to the 
simple fact that on average there may be more sources in one direction (towards the Galactic 
center) than in the opposite direction. The first moment, 5a2, as calculated above refers 
to the fluctuations on top of this non-vanishing mean anisotropy and whether it dominates 
or not again depends on the specific distribution of sources. For the source density profiles 
that will be used in the sections below, though the energy dependence of the amplitude of 
anisotropy is approximately the same for both terms, the scaling with the halo size H is 
different. The anisotropy arising from the overall inhomogeneous source distribution has a 
faster dependence on H. For small values of H, comparable with the large scale gradient 
in the source distribution, Sai is comparable in magnitude to 5a2, but becomes increasingly 
dominant with increasing H. Although the energy dependence of the mean anisotropy over 
many source realizations for the two terms is similar for our assumed T m i n (both proportional 
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to D(E)), when Sa2 is negligible compared to Sai the anisotropy amplitude approaches 
a monotonic function of energy. At the same time, the phase of the anisotropy is very 
different: it correlates strongly with the direction of the Galactic center for large H while it 
is a stochastic variable for small H, a symptom of the dominance of nearby recent SNRs in 
this latter case. The importance of nearby sources in terms of anisotropy was also discussed 
in [17] where the authors adopt a stochastic generation method for SNRs in the Galaxy 
similar to the one used in this paper. 

4 CR anisotropy for realistic distributions of SNRs 

To first approximation, the anisotropy of CRs that is observed at Earth results from the 
combination of two main physical ingredients: the global spatial distribution of SNRs in the 
Galaxy, which in turn leads to gradients in the CR distribution, and the proximity effect 
of nearby and/or recent supernova events. In general other elements can also affect the 
anisotropy, like, for instance, longer residence times of CRs close to their sources resulting 
from the self-generation of turbulence. 

In this section we summarize the two models of spatial distribution of SNRs, already 
discussed in Paper I, the cylindrical model and the spiral model. The cylindrical model is 
basically that discussed, for instance, in [18]. The other model (hereafter "spiral model") is 
instead an attempt at taking into account the spiral distribution of sources in the Galaxy: 
for this we adopt the formalism of [19]. 

In both cases the sources are assumed to have a mean radial distribution in the Galaxy 
following the function [18]: 

where (3 = 3.53 for R & = 8.5 kpc. 

The constant A is determined from the normalization condition: 

dr 2irr f(r) = 1 -> A = ^ • (4.2) 

The SNR distribution in the z direction is assumed to be 

(4.3) 

where A z = 1 is again derived from the normalization condition. 

In the cylindrical model, the positions of the sources are chosen by drawing at random 
values of r and z from the distributions above. The x and y coordinates are chosen by 
generating a random angle < eft < 2tt and using the given value of r. Supernovae are 
generated at a rate 1Z = 1/(30 — 100) yr _1 , and the generation of new sources is continued 
until a time span much larger than H 2 / 'D(E) is covered, in order to make sure that the 
stationary solution has been reached at the lowest rigidities of interest for us (~ 1 TV). 
Typically in our calculations we choose a time span of 80 million years, more than sufficient 
to ensure stationarity at rigidities larger than 1 TV. 

In the spiral model, the procedure we adopt is similar to that introduced in [19] and 
summarized in Paper I: the generation of the position in the z direction is the same as above, 
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Figure 1. A face-on view of the spatial distribution of SNRs in the Galaxy for two different models 
of the spiral structure: w — 5 kpc on the left and w = 2.8 kpc on the right. About 3 x 10 4 
are shown in each panel. Units are in kpc and the position of the Sun is at (x, y) = (Rq,0). 



therefore we will not discuss it any further. A radial coordinate r is drawn at random from 
the average distribution in Eq. 4.1; at this point a random natural number between 1 and 
4 is chosen from a flat distribution. This number identifies the arm in which the supernova 
is localized (Norma, Carina-Sagittarius, Perseus, Crux-Scutum, as in Table 1 of Paper I). 
At this point an angular position along the arm is associated to the SNR, according to the 
relation: 

9(r) = K log (j^J+e . (4.4) 

The parameters K, r$ and 6q are reported in Table 1 of Paper I. 

Following the prescription of [19], we blur the angle 6(r) by 6 CO rr exp(— f [kpc] jw [kpc]): 
here 9 corr is chosen from a flat random distribution between and 2ir, while for w we consider 
two reference values, w = 2.8 and w = 5 kpc and study how the results are affected by this 
choice. The distribution proposed in Ref. [19] has w = 2.8 kpc, but the choice of the value 
of this parameter does not appear to be strongly data motivated. 

Similarly to the angular position, the radial coordinate is also blurred by choosing a 
new value of the radius from a normal random distribution with mean f and variance 0.07f 
[19]. 

In Fig. 1 we show how the spiral structure changes depending on the value of w: the 
figure shows the distribution of ~ 30, 000 SNRs in a case where w = 5 kpc (left) and in a 
case with w = 2.8 kpc (corresponding to the value adopted by [19]). The position of the Sun 
is at (x,y) = (i? ,O). 

For each realization of source distributions in the two models the contribution to 
anisotropy is calculated in all directions, so that both a magnitude and a phase can be 
evaluated. 

4.1 Anisotropy for the cylindrical model 

In this section we discuss the anisotropy signal for the case of cylindrical symmetry in the 
SNR distribution. In Fig. 2 we illustrate the anisotropy amplitude for 10 realizations of the 
source distribution in the cylindrical model, using 5 = 1/3, H = A kpc and a rate of supernova 
explosions in the Galaxy 1Z = 1/100 yr _1 (left) and 1Z = 1/30 yr _1 (right). In all cases we 
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Figure 2. Anisotropy amplitude for ten random realizations of sources in the cylindrical model, 
assuming S = 1/3 and a SN rate TZ = 1/100 yr" 1 (K = 1/30 yr" 1 ) on the left (right). The halo 
size is H = 4 kpc. The injection spectrum is assumed to have slope (below the cutoff) such that 
7 + S = 2.67. The data points arc from [20-22]. 



impose that the slope 7 of the injection spectrum is related to 5 through 7 + 5 = 2.67, in 
order to ensure a good fit to the CR spectrum at Earth (see Paper I). The red, staircase line 
represents the average amplitude calculated using the 10 random realizations. 

In all figures the (black) crosses, the (blue) diamonds and the (orange) stars are taken 
from Ref. [20]. The (green) triangles are from EASTOP [21, 23] and the (red) squares are 
the Akeno data points [22]. The oblique (red) lines at high energy show the upper limits on 
the amplitude of anisotropy from KASCADE and GRANDE [24]. 

The comparison between the two panels shows that the spread in the anisotropy patterns 
is not affected in a significant way by the SN rate. This can be qualitatively understood if 
one considers that for H = 4 kpc, the anisotropy signal is already dominated by 5ai (see 
§ 5). Looking at Eq. 3.5 one sees that the rate of Supernova explosions TZ only enters (Jcr) 
(and the same is true for ucr) through the normalization of the probability distribution. It 
is then clear that any dependence on TZ will disappear when 6ai is obtained as the ratio 
between (Jcr) and ucr- Both panels of Fig. 2 show very clearly the strong dependence 
of the strength of anisotropy on the specific realization of source distribution, thereby also 
disproving the naive expectation that the anisotropy should be a growing function of energy 
with the same slope as the diffusion coefficient D{E). Whenever the small scale contribution 
is not negligible, the observed anisotropy can in fact even be a non monotonic function of 
energy, with dips and bumps, and with wide energy regions in which it is flat with energy, 
quite like what the data show at energies E < 10 5 GeV. It is interesting however that none 
of our realizations of the source distribution leads to anisotropies as low as the one suggested 
by the data in the energy region 10 5 - 10 6 GeV (contributed by the EASTOP experiment). 

Data in this region are in fact somewhat puzzling because they are so low as to suggest 
that the Compton-Getting effect [25] leads to a level of anisotropy close to the lowest expected 
limit. The Compton-Getting anisotropy is estimated to be between 3 x 10~ 4 and 10~ 3 
depending on the velocity with which the Earth moves with respect to the rest-frame of the 
CR scattering centers. This velocity is not known and the above estimates refer to a velocity 
range from a minimum of ~ 20 km/s to a maximum of ~ 250 km/s, corresponding to the 
motion of the solar system through the Galaxy [26] . It is clear that the measured anisotropy 
between 10 5 and 10 6 GeV is only marginally consistent with a velocity of few tens of km/s 
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Figure 3. Anisotropy amplitude for ten random realizations of sources in the cylindrical model, 
assuming S = 0.6 and a SN rate 1Z = 1/30 yr _1 . The halo size is H = 4 kpc. The injection spectrum 
is assumed to have slope (below the cutoff) such that 7 + 5 = 2.67. 

at most. 

We also checked the effects of decreasing further the source rate, which could be the 
case if the bulk of CRs does not come from standard SNe but rather from rarer events, like 
for example an especially energetic sub-sample of SNe or GRBs. The resulting anisotropy is 
somewhat larger at low energies, on average: the data can still be easily reproduced at the 
low and high energies, but the central, more problematic region is now more extended, in 
general, to the left than in Fig. 2, approximately ranging from few xlO 4 to 10 6 . 

In Fig. 2 we adopted a diffusion coefficient scaling with E 1 ^ 3 . The energy dependence 
of the diffusion coefficient is however the subject of an ongoing debate: given D(E) oc E 5 it 
is controversial whether 5 is 1/3, 1/2, 0.6 or even larger (see [27] and references therein). 

The all-particle spectrum alone, while giving some indications that 5 = 1/3 could be 
preferable (see Paper I), does not allow one to really clinch the question. This is because 
the all-particle spectrum only depends on the combination 5 + 7. In principle the B/C ratio 
would allow a direct measurement of 5, if this ratio could be measured at sufficiently high 
energies. Unfortunately at the present time the error bars on this quantity are still large 
enough to allow for ambiguity in the best fit value (see for instance [28]). 

Since the anisotropy 5a is defined as the ratio between the density gradient and the 
density, 7 does not appear in 5a while 5 does (see also expressions 3.11 and 3.13 for the 
simplified case of a uniform distribution of the sources). In Fig. 3 we plot the amplitude 
of the anisotropy computed for ten different realizations of the source distribution in the 
cylindrical model: a slope of the diffusion coefficient 5 = 0.6 is assumed, while the other 
parameters are all the same as for the plot in the right panel of Fig. 2. 

As well as in the case 5 = 1/3, also for 5 = 0.6 the amplitude of the anisotropy is a 
complex function of energy as a result of the cosmic rays contributed by nearby recent SNRs. 
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However, for 6 = 0.6 the amplitude of the anisotropy appears to be systematically larger 
than the observed one at all energies. In other words, fast diffusion leads to exceedingly 
large anisotropy which seriously challenges the models that require large values of 5 (see for 
instance the discussion in Ref. [27] for the cases in which a convective wind is included). It 
is worth noticing that at very high energies the amplitude may exceed unity. These cases 
clearly suggest that the diffusive paradigm may break down for very nearby sources of CRs, 
as already discussed in Paper I. 

We think that the results just showed provide clear evidence in favor of a diffusion 
coefficient with a weak dependence on energy. This finding is of crucial importance in several 
respects. The fact that the data suggest a value 5 = 1/3 is comforting in some respects and 
puzzling in some others, in relation to our understanding of CR acceleration and propagation. 
On the one hand, 6 = 1/3 gives the exact energy dependence of D{E) that Kolmogorov-type 
turbulence would provide, so propagation follows a framework that was not unpredicted from 
the theoretical point of view. On the other hand, however, as we already mentioned in Paper 
I, 5 = 1/3 implies that the injection spectrum should be a power law with slope 7 ~ 2.3 — 2.4 
which is much steeper than what acceleration theory has been predicting for the last several 
years, especially due to the non linear effects which make the spectrum concave. We defer 
further discussion of this point to § 6, while below we discuss the dependence of anisotropy 
on other unknowns of the problem, like the real source distribution, the size of the halo and 
the detailed time-dependence of the injection of CRs by each source. 

The results discussed so far are obtained within the assumption that each source releases 
accelerated particles in a time that is short compared to all the relevant time-scales of the 
problem; moreover, the released particles are assumed to have a power-law spectrum. As 
already anticipated it is wise to also consider a different scenario, involving time-dependent 
injection (see § 2). In this model, the release of CRs in the ISM occurs continuously during 
the Sedov- Taylor phase of expansion of the SNR, with higher energy particles being released 
earlier (close to the beginning of the Sedov phase) and particles of gradually decreasing energy 
at correspondingly later times. The adopted dependence on time and energy of the injected 
particle distribution is that illustrated in § 2, Eq. 2.6 and Eq. 2.7. It is however necessary to 
keep in mind that this recipe should not be considered as a realistic attempt to model the 
complex escape of CRs from a SNRs, an objective that is way beyond the purpose of the 
present paper and that has been discussed at length elsewhere (see e.g. [15]). Our purpose 
here is only to explore the effect of gradual injection of CRs into the ISM on the amplitude 
of anisotropy. 

In Fig. 4 we plot the anisotropy resulting from ten different configurations for two 
different choices of the SNR lifetime, tsnr = 3 x 10 4 on the left and tsnr = 3 x 10 5 on the 
right. The resulting curves are qualitatively similar to those obtained for bursting injection 
of CRs into the ISM, with possibly some additional bumpiness caused by the fact that some 
nearby sources did not have time yet to generate the low energy particles. It remains true 
that there are several realizations of source distributions which appear to provide a reasonable 
fit to the observed anisotropy amplitude. The data points between 10 5 and 10 6 GeV remain 
somewhat peculiar. 

The two different durations of the Sedov- Taylor phase do not make much difference, 
once the parameters are adjusted in such a way as to guarantee that the final slope of the 
observed spectrum, j b s , is the same: this requires that once a is determined from the 
condition a = \og(E M / 'E ) / ' log(rsNR/Ts T ), P is given as (3 = a{^ ohs - 2-5). 
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Figure 4. Anisotropy amplitude for ten random realizations of sources in the cylindrical model in 
the case of continuous escape of particles from the upstream of the accelerating blast wave. The two 
different panels are relative to different assumptions on the SNR lifetime: tsnr = 3 x 10 4 yr on the 
left; tsnr — 3 x 10 5 yr on the right. Other parameters are: S = 1/3, a SN rate 1Z = 1/30 yr -1 and 
a halo size H = 4 kpc. The injection spectrum is such as to ensure 7 + S = 2.67, see text for details. 



4.2 CR anisotropy in the spiral model 

In this section we discuss how the anisotropy of CRs is affected by the spiral structure of the 
Galaxy. 

In Fig. 5 we plot the anisotropy computed for a number of different realizations of the 
sources' distribution in the spiral model. The two panels relate to two different choices for 
the parameter w introduced above, the same two values for which the spatial distribution of 
the sources is plotted in Fig. 1, w = 2.8 kpc and 5 kpc. First thing we notice in Fig. 5 is that 
in the left panel {w = 5 kpc) the dispersion among the different curves is much larger than on 
the right {w = 2.8 kpc). The reason for the different shape of the curves for the anisotropy 
amplitude as a function of energy for the two values of w is related to the large scale structure 
of the spiral: as we discussed several times, the CR density gradient is determined by both 
the large scale distribution of SNRs and by the accidental proximity of nearby recent sources. 
Whether the anisotropy is dominated by one or the other effect depends upon whether the 
halo size H is larger or smaller than the spatial scale on which gradients in the source 
distribution appear. In the case of a spiral with a tight spread (w = 2.8 kpc), the large scale 
distribution of the sources is such that the nearby recent SNRs are subdominant and in fact 
the curves showing the amplitude of anisotropy are very regular and monotonically increasing 
with energy In the case w = 5 kpc, the sources are more spread out around the centroid 
of the arms and it becomes more likely to have recent SNRs located in the proximity of 
the Sun's location, therefore the anisotropy shows the bumps and dips that we have already 
discussed in the previous section. 

This interpretation of the results of our calculations is also confirmed by the investigation 
of the phase of the anisotropy vector. In Fig. 6 we plot &rctan(5 y /5 x ) (phase of the anisotropy 
vector in the xy plane of the Galactic disc) as a function of energy The left (right) panel 
refers to w = 5 kpc (w = 2.8 kpc). As expected, in the case of narrow arms (w = 2.8 
kpc) the anisotropy is dominated by the large scale distribution of sources and the phase has 
little spread around zero. On the other hand, for broader arms (w = 5 kpc) the fluctuations 
induced by recent nearby SNRs dominate the anisotropy amplitude and the phase has a wide 
spread depending on the realization of the source distribution. 
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Figure 5. Anisotropy amplitude for ten random realizations of sources in the spiral model. The 
curves in the two panels are obtained for different values of w, with w = 5 kpc on the left and w = 2.8 
kpc on the right (source distributions as in the left and right panel of Fig. 1 respectively). Bursting 
injection is considered. Other parameters arc: S = 1/3, a SN rate 1Z = 1/30 yr _1 and a halo size 
H = 4 kpc. The injection spectrum is such as to ensure 7 + S = 2.67, see text for details. 
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Figure 6. Phase of the anisotropy for the same cases as in the left and right panels of Fig. 5. 



It is worth stressing that the phase of the anisotropy as derived by our calculations 
shows instances of abrupt changes in contiguous energy bins, as a result of the dominance, 
in those bins, of two (or more) different sources (the same behaviour is found in the case 
of cylindrical distribution of sources). This type of behaviour is reminiscent of that found 
by [23] in the two energy bins E = 1.1 x 10 14 eV and 3.7 x 10 14 eV. In the same paper the 
authors also tried to make sense of the slope of the amplitude of the anisotropy between the 
same energy bins, and found a value 5 ~ 0.74, which they interpreted in terms of energy 
dependence of the CR diffusion coefficient. As showed by our calculations, such a slope is 
unlikely to have anything to do with the diffusion coefficient, and is rather the result of the 
large fluctuations induced by local sources. 



5 The dependence of the anisotropy amplitude on the size of the halo H 

As we discussed in § 3, for a homogeneous distribution of sources in a thin disc the only 
contribution to the anisotropy is the one due to the proximity of local recent sources. In 
fact, the non central position of the Sun in the disc also gives rise to a small anisotropy but 
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Figure 7. Relative importance of the large scale (Sai, solid lines) and small-scale (6a2, dashed lines) 
contribution to the anisotropy signal. The two terms are plotted for different sizes of the halo (H = 1 
kpc (red), H = 2 kpc (purple), H = 10 kpc (blue). The diffusion coefficient is assumed to scale as 
E 1 / 3 . 



this component is negligible compared to that associated with fluctuations. Moreover, the 
anisotropy as derived in § 3 (see Eq. 3.13) is independent of the size H of the halo, once the 
constraint D2$/Hk pc ~ 1 is taken into account. It follows that making the size of the halo 
bigger should not lead to any appreciable change in the anisotropy amplitude. 

This situation changes however when the realistic distribution of sources in the disc is 
taken into account. 

We computed the terms 6ai and $A2 according to Eqs. 3.3 and 3.4 where in the definition 
°f (Jcr) and {SJcrSJcr} (Eqs. 3.5 and 3.6) we used the probability distribution described 
by Eq. 4.1 and normalized to the supernova rate: in practice we computed the anisotropy 
corresponding to a continuous (instead of stochastic) version of the cylindrical model for 
source distribution. The results of this calculation are plotted in Fig. 7. 

It is clear from this figure that the scaling with H of the anisotropy as due to the large 
scale distribution of sources is different from that due to fluctuations. In particular the former 
contribution increases appreciably with increasing halo size, while the small-scale term is left 
more or less unaltered: as a result the small scale contribution, while comparable to 5ai 
for values of H that are smaller than the scale of the gradient in the source distribution, 
Rq/(3 ~ 2.5 kpc (see Eq. 4.1), becomes close to negligible for H = 10 kpc. 

In order to illustrate this effect within our model of stochastic source distribution we 
recalculate the amplitude of the anisotropy for 5 = 1/3 and 1Z = 1/100 yr" 1 using a smaller 
halo size with H = 1 and H = 2 kpc. The results are illustrated in Fig. 8: both panels 
refer to the cylindrical model, with H = 1 kpc on the left and H = 2 kpc on the right. One 
can easily realize that increasing the size of the halo reduces the fluctuations. This can be 
understood because, as we stressed above, for large values of H the anisotropy amplitude 
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Figure 8. Anisotropy amplitude for ten random realizations of sources in the cylindrical model for 
H = 1 kpc (left panel) and H = 2 kpc (right panel). 




Figure 9. Phase of the anisotropy in a cylindrical model with H = 1 kpc (left panel) and H = 4 kpc 
(right panel). 



is dominated by the regular, large scale, distribution of the sources. In both cases, most 
realizations lead to an anisotropy that mostly grows with energy. In other words, for large 
values of H it becomes harder to explain the anisotropy signal observed at Earth. In this 
sense, though not ruled out, large values of H appear to be disfavored. Interestingly enough, 
in the case with H = 1 kpc, there are a few realizations that lead to an anisotropy which is 
remarkably similar to the observed one. 

The trend just described can be illustrated more clearly by using the phase of the 
anisotropy vector, as plotted in Fig. 9 for a halo of size H = 1 kpc (left) and H = 4 kpc 
(right). For small values of H the phase varies wildly reflecting the occasional dominance of 
a nearby recent source. Again, this behavior is reminiscent of that found by the EASTOP 
experiment [23], as discussed above. For H = 4 kpc the main contribution to anisotropy 
comes from the inhomogeneous source distribution in the Galactic disc, and the energy 
dependence of the phase of the anisotropy becomes much more regular, with an offset with 
respect to zero that reflects the presence of some nearby source. 
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6 Conclusions 



We determined the amplitude of the anisotropy expected in the CR flux at Earth if SNRs 
are their main sources. We considered different scenarios in terms of the spatial distribution 
of Galactic SNRs, of the injection of CRs into the ISM from individual SNRs, and of the 
energy dependence of the diffusion coefficient D(E) cx E s . 

The amplitude of the anisotropy reflects the gradient in the CR density at Earth and 
includes both a contribution from the large scale distribution of sources in the Galaxy, and 
from the random nature (in space and time) of nearby recent SNRs. Whether the first or 
the second contribution dominates depends on the specific source distribution in space and 
on the choice of the halo size. For values of H < 4 kpc the predicted anisotropy is the 
result of comparable contributions from the fluctuations induced in the CR flux from nearby 
sources and the inhomogeneous source distribution. This implies that the curve describing 
the anisotropy as a function of energy is characterized by bumps and dips, and is, in general, 
all but monotonic. In other words the naive and often quoted expectation that the anisotropy 
scales as cx D(E) is far from what one can measure for any specific realization of the source 
distribution. When the average of the anisotropy amplitude is computed over many source 
realizations, this turns out to be closer to the expected scaling, but the agreement is still not 
perfectly realized because of the presence of nuclei. These, after propagation, have spectra 
that are different, in general, from those of protons (see Paper I). As a result, the scaling of 
anisotropy with energy, as averaged over many source realizations, turns out to be somewhat 
slower than oc D(E). 

The bumps and dips mentioned above often resemble plateau regions in limited energy 
intervals, quite similar to what is actually observed at the Earth. Although it is not possible 
to explain the detailed structure of the anisotropy amplitude as reported in the figures above, 
since this depends upon the specific realization of nearby recent SNRs in which we live, we 
claim that the qualitative features of the observed data points can be understood as resulting 
from the stochastic nature of SNRs as sources of cosmic ray protons and nuclei. 

The most important conclusion of our calculations is that the anisotropy is a strong 
function of the diffusion coefficient. Assuming that the ratio H/D(E) is normalized in such 
a way as to fit the observed ratio of B/C at ~ 10 GeV/n, we showed that the amplitude of 
the anisotropy for 5 = 0.6 exceeds the observed amplitude by up to two orders of magnitude, 
depending on the realization of sources and on the energy. Larger values of 5 make the 
problem even more severe. 

This is an important conclusion in that 5 > 0.6 would appear to be required by both the 
standard (test-particle) theory of particle acceleration at strong shocks (in which the slope 
of the injected spectrum is ~ 2) and by the non-linear version of the theory, in which at 
high energy it is usually found that the injection spectrum is even harder than E~ 2 (see for 
instance [29]). Only very recently [6, 8, 9] it has been proposed that steeper spectra might 
result from phenomenological (and in general not very well justified in terms of basic physical 
principles) approaches to the problem of determining the velocity of the scattering centers. 
The fact that we find large values of 5 to be incompatible with the measured amplitude of 
the anisotropy is an additional argument in favor of the existence of physical phenomena that 
steepen the injection spectrum, with all the implications that this finding involves. 

All the conclusions relative to the choice of the diffusion coefficient remain valid also 
when the spiral structure of the Galaxy is taken into account. However, the spiral structure 
in the spatial distribution of SNRs affects the contribution to the amplitude of the anisotropy 
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as due to the large scale gradient of CRs. The scale of the spiral arms is found to have an 
important effect on the anisotropy. If the distribution of the sources around the centroid 
of the arms is very narrow (case w = 2.8 kpc above), the gradient induced by the spiral 
structure on the CR density dominates upon the fluctuations, hence the amplitude of the 
anisotropy grows with energy rather regularly, roughly as (but slightly more slowly than) 
oc D{E). For broader arms (w = 5 kpc above), the role of fluctuations becomes dominant 
again, so that the amplitude shows the standard bumpy structure discussed earlier. 

In both cases of cylindrical and spiral distribution of SNRs in the Galaxy, the pattern 
of anisotropy is weakly affected by the size of the halo, H. However, when H becomes much 
larger than the scale of gradients, namely Rq/(3 in the cylindrical model, and w in the spiral 
case, the contribution due to the CR gradient induced by the large scale source distribution 
dominates upon the fluctuating terms. Increasing the value of H in general leads to a very 
regular trend of nearly monotonic growth of the anisotropy amplitude with energy, quite 
unlike what is observed. Such a behavior also reflects on the phase of the anisotropy which 
is very irregular for small H and very regular for larger haloes. This argument plays against 
models with large halos, which however appear to be favored in calculations carried out with 
GALPROR 

We conclude our discussion of the results of our calculations with a short comment 
on the importance of using a stochastic method like the one discussed here in order to 
have a reliable determination of the amplitude (and phase) of the anisotropy, rather than a 
propagation code such as GALPROP or DRAGON. These latter approaches can only pick the 
anisotropy (proportional to the CR gradient) induced by the large scale distribution of SNRs 
(or other sources) in the disc of the Galaxy, while missing completely the fluctuations, that 
however represent the dominant contribution to anisotropy for most values of the relevant 
parameters. At the same time it is fair to recall that both numerical propagation codes as 
well as the calculations presented here are not able to treat a major complication intrinsic 
to the problem, namely the fact that the propagation in the Galactic magnetic field might 
be more complex than modeled through a spatially constant, isotropic diffusion coefficient. 
For instance, it is very likely that the amplitude and phase of anisotropy are affected by 
anisotropic diffusion, namely preferential diffusion parallel to the magnetic field lines. On 
the other hand it is plausible that this effect is somewhat mitigated by the random walk of 
magnetic field lines, which leads to make diffusion close to isotropic. This also should reflect 
in a more erratic behaviour of the phase of the anisotropy vector, which seems pretty similar 
to what the data show. 
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